Analytische oefeningen t.b.v de Voortoets

@Theo Olsthoorn (12-03-2024 .. 30-11-2024)

Intro, verantwoording

Na de vorige vergadering van 12 maart 2024 heb ik een aantal aspecten geanalyseerd die van belang zouden kunnen zijn voor een relatief eenvoudige, analytische analyse van de impact van een ingreep op een grondwatersysteem. De analyse heb ik uitgevoerd en gedocumenteerd in voorliggend Jupyter (Python) notebook, dat uitleg en code bevat om het uitgelegde te kwantificeren en te laten zien.

Voor een aantal complexere concepten zoals vertraagde nalevering, debietverloop van een bemaling met constante verlaging, hoe lang het duurt voor de verlaging stationair wordt e.d. zijn vereenvoudigingen voorgesteld die in de praktijk goed zullen werken.

Voor een zinvolle analyse van een fysische ingreep in het grondwatersysteem, zoals een nieuwe puntonttrekking, bemaling of verandering van de loop of het niveau van oppervlaktewater, is een beeld nodig van de opbouw van de ondergrond en van de drainage in het gebied (de randvoorwaarden). Bovendien is men geïnteresseerd in de overlap van de impact met bijzondere beschermingsgebieden. Deze drie vormen van informatie zijn gebonden aan beschikbare kaarten zoals die van de habitatgebieden, het oppervlaktewater, en bodemlagen. De laatste twee zijn aanwezig in de kaarten waarop het grondwatermodel van Vlaanderen is gebaseerd. Deze kaarten bevatten ook de benodigde informatie over de waarden van de bodemconstanten, zoals doorlaatvermogen, weerstand tussen lagen en bergingscoëfficiënten. Mogelijk kunnen op basis van de beschikbare laagverbreidingskaarten ook complexere situaties worden gesignaleerd, zoals breuken en andere scherpere overgangen tussen gebieden, die niet met eenvoudige formules kunnen worden geanalyseerd of waarvoor een aangepaste berekening kan worden voorgesteld of voorgeschreven. Het uit kaarten halen van de randvoorwaarden voor een berekening is vermoedelijk het meest complex of valt in de praktijk het meest op af te dingen. Uiteraard kunnen randvoorwaarden worden voorgeschreven zoals een invloedsradius of duur van de onttrekking waarmee moet worden gerekend, zoals dat nu reeds in de Voortoets het geval is.

Voor de benodigde onderliggende gegevens moet er i de Voortoets toegreep zijn tot het Vlaamse grondwatermodel, of althans de kaarten waar dit op gebaseerd is, zodat deze kaarten als een ruimtelijke database kunnen worden beschouwd en bevrraagd op bodemconstanten die voor een gegeven locatie moeten worden gebruikt.

Door op een dergelijke manier de voor elke vraag benodigde informatie op te vragen, kan de onderliggende machinerie van de voortoets steeds gemakkelijk worden aangepast en verbeterd naar nieuwe inzichten.

Het resultaat van de Voortoets zol op deze wijze ook steeds in lijn zijn met de eventueel naderhand uit te voeren bredere analyse, waar dan zonodig een ruimtelijk grondwatermodel aan te pas komt, dat immers op dezelfde gegevens is gebaseerd. Een kernpunt zou dus moeten zijn dat de gegevens die gebruikt worden in de voortoets dezelfde zijn als die in het ruimtelijke model van Vlaanderen, waarbij de voortoetsberekeningen zich zullen baseren op de ondergrond gegevens ter plaatse van de ingreep en het ruimtelijk model met de ruimtelijke variatie rekening houdt.

De info-vraag van de Voortoets zal altijd zeer beperkt zijn, zodat de ICT-belasting navenant laag blijft en de informatie dus real-time over het internet (via een URL) moet kunnen worden opgevraagd en worden opgezocht op de ruimtelijke kaarten.

Imports voor de benodigde functionaliteit

Basisfuncties die verderop worden gebruikt:

Verlaging binnen een cirkelvormmige bouwput met uniforme onttrekking langs de omtrek

We gaan uit van de meest eenvoudige formule voor de verlaging, namelijk Dupuit. De onttrekking langs het scherm met radius $\rho$ is gelijk aan $dq = \frac{Q_0}{2 \pi \rho} \rho d\theta = \frac{Q_0}{2 \pi}d\theta$. De verlaging hierdoor in een willekeurig punt in de bouwput, gegeven door coordinaten $x_p=a\rho$, $y_p=0$, met $0 \le a \le 1$, waarbij het centrum van de bouwput coordinaten $0, 0$ heeft. De radius waarvoor de verlaging gelijk is aan nul is veel groter dan die van de boutput en is gelijk aan $R$.

De afstand $r$ tussen het punt op de omtrek van de bouwput en het beschouwde punt is gelijk aan

$$r = \rho\sqrt{\sin^2\theta + (\cos\theta)^2 - a}=\rho\sqrt{1+a^2-2a\cos\theta}$$

De verlaging $ds$ door onttrekking van het stukje scherm gegeven door $d\theta$ is gelijk aan

$$s = \frac {1}{2 \pi}\frac{Q_0}{2 \pi kD}\intop_0^{2 \pi}\ln[\frac \rho R \sqrt{1 + a^2 - 2 a \cos\theta}] d \theta$$$$s = \frac{Q_0}{2 \pi kD}\ln \frac r R + \frac {1}{2 \pi}\frac{Q_0}{2 \pi kD}\intop_0^{2 \pi}\ln[\sqrt{1 + a^2 - 2 a \cos\theta}] d \theta$$

De tweede term is lastig analytisch te integreren, maar bij numerieke integratie blijkt deze altijd gelijk aan nul te zijn. Met andere woorden, de verlaging binnen een cirkelvormige bouwput die unirfoem vanaf de rand wordt bemalen is altijd uniform en gelijk aan de verlagin op de rand van de bouwput veroorzaakt door een onttrekking in het centrum ervan.

Hieronder laten we zien dat de integraal nul is en daarmee ook de tweede term.

Dus ook in de situatie met lek, De Glee, is dit het geval zolang de radius van de bouwput $\rho << \lambda=\sqrt{kD c}$, de spreidingslengte van de semi-gespannen aquifer. Immers dan geldt

$$s \approx \frac{Q}{2 \pi kD} \ln \frac{1.123 \lambda}{\rho} = \frac{Q}{2 \pi kD}\ln \frac{R}{\rho},\,\,\,\,\,met\,\,\,\,R=1.123 \lambda$$

In de niet-stationaire situatie volgens Theis, voldoen de stijghoogteverschillen als functie van de afstand aan die van de formule van Dupuit. Dit impliceert dat, na voldoende tijd, dat wil zeggen wanner de invloedsstraal ruim groter is dan de straal van de bouwput, ook dan geldt dat de verlaging binnen de cirkelvormige bouwput exact vlak. Deze verlaging is dan bovendien gelijk aan de verlaging op de rand van een bouwput door een even grote ontrekking in het centrum van de bouwput. Kort na de start van de bemaling is nog niet het geval.

Het voorbeeld hieronder demonstreert dit.

De conclusie is tweeledig.

  1. De verlaging binnen en cirkelvormige bouwput is overal gelijk bij een schermonttrekking langs de omtrek. En dus bijna overal gelijk bij onttrekking middels een aantal putten verdeeld langs de omtrek.
  2. De verlaging binnen de bouwput bij onttrekking langs de omtrek is gelijk aan de verlaging op de omtrek bij onttrekking met een put met hetzelfde totale debiet in het centrum van de bouwput.
  3. Dit laatste is exact het geval in de stationaire situatie, doch de fount in de niet stationaire situatie is al snel gering, in feite te verwaarlozen wanneer de tijd voortschrijdt.

Hierna wordt aangetoond dat de verlaging binnen de bouwput (na enige tijd) vlak is en gelijk aan de verlaging op de rand van de bouwput door een put in het centrum van de bouwput met hetzelfde totale debiet.

Wanneer het scherm geen onttrkkingsscherm is, maar de bemaling wordt gedaan met een aantal putten op een circle rond het hard van de bouwput dan geldt

$$s = \frac{Q_0}{ 2\pi kD} \ln \frac R r$$

N bij N putten alle op de cirkel met straal $R$ geldt voor het punt in het midden van de cikel

$$s(r=0) = \frac{Q_0}{ 2\pi kD} \ln \frac R r = \frac{\sum_{i=1}^NQ_i}{2 \pi kD} \ln \frac R r$$

Dus bij bemaling met N putten langs de omtrek van een cikelvormige bouput is de verlaging in het centrum van deze cirkel gelijk aan die op afstand $R$ veroorzaakt door het pompen in het centrum van de cirkel met een debiet gelijk aan de som van de debieten van de bronnen op de cirkel. Het maak niet uit waar de putten op de cirkel staan en hoe groot hun onderlinge debieten zijn.

Het voorbeeld hierna geeft de verlaging op langs rand van de cirkelvorminge bouwput bij onttrekking met een (willekeurig) aantal putten, dti voor verschillende momenten in de tijd. Alleen wanneer de onttrekking uniform wordt verdeeld over de gehele cirkel valt de verlaging door de scheronttrekking op de cirkel samen met die door een enkele put in het centrum van de bouwput. In het voorbeeld kan het aantal putten worden gekozen waarover de totale ontterkking wordt verdeeld over de rand van de boutput.

De nalijlende verlaging

Hoeveel blijft er na verloop van tijd nog over van een tijdeijke onttrekking in het verleden? Deze vraag is eigenlijk alleen relevant voor de Theis situatie met freatische berging, omdat de semi-gespannen situatie volgens Hantush al snel stationair wordt en dan ook snel geen restantvergelaging meer heeft nadat de onttrekking is gestopt.

De verlaging door een onttrekking in een situatie zonder randvoorwaarden wordt goed beschreven door Theis

$$ \frac Q {2 \pi kD} W(u), \,\,\,\, u =\frac{r^2 S }{4 kD t} $$

Er is in een dergelijke Theis situatie na beëndiging van de onttrekking nog zeer lange tijd sprake van een resterende verlaging, die samen met andere niet permanente onttrekkingen toch tot een cumulatief effect zal leiden.

De resterende verlaging is dan gelijk aan $$ s(r, t) = \frac{Q}{4 \pi kD} \left[ W_t \left(\frac{r^2 S}{4 kD t}\right) - W_t \left( \frac{r^2 S}{4 kD (t - \Delta t)} \right)\right]$$

Het voorbeeld hieronder laat de verlaing als functie van de tijd zien, inclusief het naijlen ervan nadat de bemaling is gestopt.

Het cumulatieve effect van gelijktijdige en in de tijd verschoven onttrekkingen

Bij een groot aantal bemalingen willekeurig in de tijd ontstaat zo een cumulatief effect.

Veronderstel een stad met een oppervlak van 3x3 km, waar op eerste van elke maand op een willekeurige locatie in de stad een bemaling start van 1800 m3/d start, die drie maanden aanhoudt. Dat betekent dat vanaf maand 3 er steeds 3 bemalignen actief zijn elke op een willekeurige locatie in de stad. De bemalingen vinden plaats over een periode van 120 maanden. De periode daarna laat naijling zien. De periode gedurende de beschouwde 120 maanden laat het cumulatieve effect van de bemalingen zien.

De bodemconstantn $kD$ en $S$ zijn zodanig dat bij 1800 m3/d op 10 m afstand na 120 dagfen een verlaging optreedt van 3 m. Dit zijn allemaal redelijke aannamen zoals je die in een stad zou verwachten.

Voorbeeld cumulatief effect van over een stad verdeelde gelijktijdige en in de tijd verschoven bemalingen

Deze denkbeeldige stad van 3x3 km heeft 100 bemalingen verdeeld over 10 jaar, die alle gedurende 3 maanden actief zijn. De locaties zijn willekeurig over de stad verdeeld. Elke maand start een nieuw onttrekking op een van de 100 locaties, en duurt da 3 maanden. Dit impliceert dat er vanaf maand 3 steeds 3 bemalingen actief zijn op 3 van de 100 willekeurige locaties. Na 103 maanden zijn alle bemalingen beëindigd.

Locaties van de bemalingen

Het cumulatieve effect van de bemalingen

Formule van Verruijt versus die van Dupuit en de formule van Blom

Het idee achter deze formules

Het idee achter de formule van Verruijt is een onttrekking in het centrum van een circulair gebied dat gevoed wordt met een constant en uniform neerslagoverschot. Het stroombeeld is hiermee uniek bepaald. Het stijghoogtebeeld is echter alleen bepaald wanneer deze ergens, op een, op zichzelf willekeurige afstand van de put wordt gefixeeerd. De formule van Verruijt is stationair and kan direct worden afgeleid voor een watervoerend pakket met vrije watertafel, waarin de effectieve dikte van het watervoerend pakket varieert met de hoogte van de watertafel. Als het gebied groot genoeg is, is er altijd een afstand waarbinnen de voeding van de put geheel richting put stroomt, en waarbuiten dit van de put af stroomt. Deze afstand vormt de waterscheiding, waarop de watertafel horizontaal is. De afstand tot deze waterscheiding neemt toe met de grootte van de onttrekking.

De wiskundige oplossing voor de situatie die Verruijt voor ogen had, is direct vergelijkbaar met die volgens Dupuit, die geen rekening hield met het neerslagoverschot. De verlaging volgens Verruijt en Dupuit, dat wil zeggen de oorspronkelijke stijghoogte minus de nieuwe, verlaagde stijghoogte, is voor beide situaties dezelfde.

De situatie die Blom voor ogen had, is een put in heen gebied met voldoende sloten, zodat de weerstand de sloten mag worden beschouwd als een vlakdekkende drainageweerstand. Voorafgaand aan de onttrekking draineren deze sloten de voeding van het gebied vanuit het neerslagoverschot. Binnen een nader te bepalen afstand rond de put vallen alle sloten droog; daarbuiten blijven zijn draineren. De drainage door de sloten die droogvallen komt volledig ten goede aan het door de put onttrokken water. Buiten het gebied met de nu droge sloten is de drainage afgenomen maar niet tot nul, en komt zodoende alleen de afgenomen drainage ten goede aan het onttrokken water. De drainage op de grens van het gebied met droge sloten is juist gelijk aan nul, daarbuiten neemt de deze asymptotisch toe tot het totale neerslagoverschot.

Bij een vlakdekkende drainageweerstand $c$, per definitie gelijk aan de gemiddelde grondwaterstand boven het slootpeil, $h$ gedeeld door het neerslagoverschot $N$, is de verlaging $s_R$ op de grens van het gebide met drooggevallen sloten en niet-drooggevallen sloten gelijk aan $s_R=Nc$.

Verruijt voor 1D stroming alleen langs de x-as

We kunnen de Formule van Verruijt zowel afleiden voor radiale situatie als die met uitsluitend stroming in de $x$-richting. Beide formules kunnen in de praktijk van pas komen en onderlinge vergelijking geeft een beter inzicht in het gedrag dat ze beschrijven.

In de eendimensionale situatie is de onttrekking geen put maar ene lijnonttrekking en is de cirkelvormige rand een rechte rand met gegeven grondwaterstand op een gegeven afstand $L$. De lijnonttrekking $Q$ heeft dimensie [L2/T] in plaats van [L2/T] bij de axiaal symmetrische situatie.

De waterbalans op afstand $x$ van de onttrekking bij neerslagoverschot $N$ is dan

$$ Q(x) = Q(0) - N x = + kh \frac{dh}{dx}$$

waarbij de stroming naar de onttrekking toe (links) posiftief wordt genomen tegen de richting $x$ in, zodat een onttrekking postief is.

$$ Q_0 x - \frac 1 2 N x^2 = \frac 1 2 k h^2 + C$$

en met gegeven stijghoogte $H$ op $x=L$ volgt voor $C$

$$ Q_0 L - \frac 1 2 N L^2 = \frac 1 2 k H^2 + C$$

Deze twee vergelijkingen van ekaar aftrekken elimineert constante $C$, zodat

$$ Q_0 (x - L) - \frac 1 2 N (x^2 - L^2) = \frac 1 2 k (h^2 - H^2) $$

oftewel

$$ h^2 - H^2 = \frac N k (L^2 - x^2) - \frac {2 Q_0} k (L - x)$$

Water divide

$$ 2 h \frac{dh}{dx} |_{x=x_D} = - 2 \frac N k x_D + 2 \frac{Q_0} k = 0. \,\,\,\,\,\,\rightarrow\,\,\,\,\,\, x_D =\frac {Q_0} N$$

Bij constant doorlaatvermogen kunnen we $h + H$ afsplitsen en schrijven als $2D$ met $D$ de pakketdikte

$$ h^2 - H^2 = (h - H)(h + H) = (h - H) 2D $$

Zodat in bij uniform doorlaatvermogen $kD$ geldt

$$ h - H = \frac N {2kD} (L^2 - x^2) - \frac {Q_0} {kD} (L - x)$$

Met $h-H = -s$ de verlaging, volgt voor deze verlaging de volgende uitdrukking

$$ s = \frac {Q_0} {kD} (L - x) -\frac N {2kD} (L^2 - x^2)$$

Blom, ééndimensionaal

Bij blom is de stijghoogte op $x=L$ niet gefixeerd, maar is de drainage daar net nul. Er is daar dus nog een restverlaging die wordt veroorzaakt door de lek (afvoerreductie) die de onttrekking veroorzaakt in het gebied waar de sloten nog wel blijven drainenen, zij het niet het volledige neerslagoverschot. Schrijven we de grondwaterstand op $x=\infty$ gelijk aan $H_\infty$ en die op $x=L$ als $H_L$ dan volgt met $h^2 - H_L^2 = (h^2 - H_\infty^2) - (H_L^2 - H_\infty^2)$

Bij variabele pakketdikte

$$ h^2 - H_\infty^2 = \frac N k (L^2 - x^2) - \frac {2 Q_0} k (L - x) + (H_R^2 - H_\infty^2)$$

De afgeleide naar $x$ is

$$ 2 h_L \frac{dh}{dx} = - 2 \frac {N L} k + \frac {2 Q_0} k$$$$ \frac{dh}{dx} = - \frac {N x} {k h} + \frac {Q_0} {k h}$$

en op $x=L$ is dit

$$ \frac{dh}{dx}|_{x=L} = - \frac {N L} {k h_L} + \frac {Q_0} {k h_L}$$

Bij vaste pakkedite $D$

Hier geldt $h^2 - H_\infty^2 = -s\,2D$ en $H_R^2 - H_\infty^2 = -s_L\,2D$

Zodat in dat geval

$$ h - H_\infty = \frac N {2kD} (L^2 - x^2) - \frac {Q_0} {kD} (L - x) + (H_R - H_\infty)$$$$ -s = \frac N {2kD} (L^2 - x^2) - \frac {Q_0} {kD} (L - x) - s_L$$

en tenslotte

$$ s = \frac {Q_0} {kD} (L - x) - \frac N {2kD} (L^2 - x^2) + s_L$$

Met als afgeleide

$$ \frac{ds}{dx} = -\frac {Q_0} {kD} + \frac {Nx} {kD}$$

en voor $x=L$

$$ \frac{ds}{dx}|_{x=L} = -\frac {Q_0} {kD} + \frac {NL} {kD}$$

Behalve het teken is de wiskundige vorm van de afgeleide hetzelfde voor de situatie met als voor die zonder vaste pakketdikte. Bij variabele pakketdikte moet in plaats van de gehele pakketdikte $D$ die ter plekke van $x=L$ moet worden gebruikt $h_L$ worden ingevuld. In veel situaties is het verschil klein. Wanneer dit verschil relevant is, moet het verhang in de situatie met variabele pakketdikte iteratief worden berekend.

Voor $x > L$, waar voeding door lek (verminderde afvoer) vanuit de sloten optreedt.

Voor $x > L$ hebben we voeding door lek vanuit de sloten via de drainageweerstand, die wordt gekarakteriseerd door de spreidingslengte $\lambda=\sqrt{kD c}$ met $c$ de drainageweerstand. De drainageweerstand is per definitie gelijk aan de gemiddelde grondwaterstand in het gebied met de sloten gedeeld doo de voeding, dus $c = (h_{gemiddeld} - H_{sloot}) / N$, met dimensie [L].

De oplossing voor de eendimensionale stationaire stroming in het gebied met nog drainerende sloten is dan

$$ s_{x>L} = s_L \exp(-\frac {x - L} \lambda),\,\,\,\,\,\, \lambda = \sqrt{kD c} $$

en de stroming

$$Q_{x \ge L} = s_L \frac{kD}{\lambda} \exp \left(-\frac {x-L} \lambda \right) $$

Onttrekking Q_0 [m2/d] op $x=0$, met voeding uit neerslag gelijk aan N, waarbij de sloten tot $x=L$ droogvallen en daarbuiten bijven draineren. Waar voor de onttrekking alle sloten de voeding wegdraineerden, is dat voor $x < L$ niet meer het geval. Deze voeding komt nu ten goede aan de onttrekking. De situatie is dan hetzelfde als in het eendimensionale geval van Verruijt.

Bij constante aangenomen doorlaatvermogen $kD$ geeft dit de volgende verlaging

$$ s = \frac{Q_0 (L - x)}{kD} - \frac{N \left(L^2 - x^2\right)}{2 kD} + s_L$$

met $s_L$ de verlaging op $x=L$

De afgeleide op $x=L$

$$ \frac{ds}{dx}|_{x=L} = -\frac{Q_0 L}{kD} + \frac{NL}{kD}$$

Voor $x > L$ wordt het pakket uitsluitend gevoed door lek, die wordt gekarakteriseerd door de spreidingslengte

$$\lambda = \sqrt{kD c}$$

waarin $c$ de zogenoemde drainageweerstand. Dit is de gemiddelde grondwaterstand tussen de sloten gedeeld door het het gedraigneeerde neerslagoverschot. De verlaging voor $x>L$ is dan

$$s_x = s_L \exp \left(- \frac {x-L} \lambda \right)$$

En het debiet $Q_{x \ge L} = -kD \frac{ds}{dx}$. Merk op dat $Q$ steeds positief is genomen in de negatieve $x$-richting, dus wanneer de stijghoogte met $x$ toeneemt en de verlaging met $x$ afneemt

$$ Q_{x \ge L} = -s_L \frac{kD} \lambda \exp \left(- \frac {x-L} \lambda \right)$$

end dus met voor $x=L$, waar $Q_L = Q_0 - L N$

$$ Q_0 - L N = -s_L \frac{kD} \lambda$$

We kennen ook de verlaging $s_L$ voor $x=L$, want daar is de voeding gelijk aan de lek via de drainageweerstand

$$x=L \,\,\,\,\,\,\rightarrow\,\,\,\,\,\, N=\frac{H_\infty - H_L}{c} = \frac {s_L}{c}\,\,\,\,\,\,\rightarrow\,\,\,\,\,\, s_L = N c$$

met $ N c kD / \lambda = N c \lambda$ volgt

$$ Q_{x \ge L} = -N \lambda \exp \left(- \frac {x-L} \lambda \right)$$

en dus hebben wel voor $x=L$

$$Q_0 - L N = N c \frac{kD} \lambda$$$$Q_0 - L N = N \lambda$$

en tenslotte

$$L = \frac{Q_0}{N} - \lambda$$

Voor de situatie met variabele dikte kan een correctie hierop worden uitgevoerd

$$L = \frac{Q_0}{N} - \frac {kh_Lc} \lambda = \frac{Q_0}{N} - \frac{k h_L c}{\sqrt{kDc}} = \frac{Q_0}{N} - \lambda \sqrt{\frac {h_L} D} $$

Voor drainageweerstand c=0, gaat de formule van Blom over in die van Verruijt

Voor $c=0$ verloopt de voeding vanuit de sloten zonder enige weerstand. De plossing is dan dezelfde als met een vaste rand op $x=L$. Met toenemende weerstand en of doorlaatvermogen neemt $L$ af. $L=0$ wanneer $\frac{Q_0} N = \lambda$ dus wanneer

$$Q_L = Q_0 = N \lambda$$

Dit is de situatie wanneer de verlaging op $x=0$ exact gelijk is aan $N c$ en de sloot op $x=0$ dus net niet meer draineert.

En wanneer $L=0$ geldt in feite dat de verlaging op $x=0$ kleiner of gelijk is aan $N c$ en de sloot op $x=0$ nog wel (enigszins) draineert.

$$ s_{x>0} = s_0 \exp \left(-\frac x \lambda \right)\,\,\,\,\,\, met\,\,\,\,\,\, Q_0 = s_0 \frac {kD} \lambda \le N c \frac{kD} \lambda = N \frac {\lambda^2} \lambda = N \lambda$$

en daar hier $s_0 = N c$ volgt

$$ Q_0 = N c \frac {kD} \lambda = N \lambda$$

Dus, $L=0$ wanneer de onttrekking zodanig is dat de verlaging op $x=0$ precies gelijk is aan de voeding.

Wellicht is interessant op te merken, dat de totale toestroming vanuit het gebied met nog wel drainerende sloten, $r \ge R$ gelijk is aan $N \lambda$ dus als het ware geschiedt vanuit een strook sloten ter breedte $\lambda$.

In de situatie met sloten die deels droogvallen kan de afstand $L$ tot waar dat het geval is direct worden berekend uit

$$L = \frac{Q_0}{N} - \lambda$$

Zoals we zullen zien is dit is bij axiaal symmetrische stroming niet het geval.

Verruijt axiaal symmetrisch

Variabele pakketdikte $h$

De formule van Verruijt voor axiale stroming naar een put gaat er ook van uit dat de onttrekking geheel wordt gevoed vanuit het neerslagoverschot. Ook hier is de stroming volledig bepaald door de onttrekking en het neerslagoverschot

$$ Q_r = Q_0 - \pi r^2 N = 2 \pi r k h \frac{dh}{dr} = \pi r k \frac {dh^2}{dr} $$

oftewel

$$ \frac{dh^2}{dr} = \frac {Q_0}{\pi k r} - \frac {r N }{k}$$

Integratie levert

$$ h^2 = \frac{Q_0}{\pi k} \ln r - \frac N {2 k} r^2 + C $$

De constante invullen geeft

$$ h_R^2 = \frac{Q_0}{\pi k} \ln R - \frac N {2 k} R^2 + C $$

Beide vergelijkingen van elkaar aftrekken elimineert deze integratieconstante

$$ h^2 - h_R^2 = -\frac{Q_0}{\pi k} \ln \frac R r + \frac N {2 k} \left(R^2 - r^2\right) $$

de afgeleide van $dh/dr$ is

$$ \frac{dh}{dr}=\frac{Q_0}{2 \pi k h} \frac 1 r - \frac{N r}{ 2 k h} $$

In deze uitdrukking komt de nog onbekende pakketikte $h$ in het rechter lid voor, maar valt eruit op de waterscheiding, $r=R_S$, wanneer we de afgeleide gelijk aan nul stellen. dit levert

$$ \frac{Q_0}{2 \pi k h} \frac 1 R_S = \frac{N R_S}{k h} $$

Oftewel

$$ R_S = \sqrt{\frac{Q_0}{\pi N}} $$

Constante pakkedikte $D$

En voor een constante pakketdikte $D$, met $h^2 - h_R^2 = (h - h_R) \, 2 D$ en $h-h_R = -s$, wordt de verlaging

$$ h - h_R = -s =-\frac{Q_0}{2 \pi kD} \ln \frac R r + \frac N {4 kD} \left(R^2 - r^2\right) $$

En dus geldt

$$ s = \frac{Q_0}{2 \pi kD} \ln \frac R r - \frac N {4 kD} \left(R^2 - r^2\right) $$

De afgeleide $ds/dr$ is

$$ \frac {ds}{dr} = -\frac{Q_0}{2 \pi kD} \frac 1 r + \frac {N r}{2 kD} $$

Het teken van de twee termen in het rechter lid zijn dan precies omgekeerd. De afgeleide gelijk aan nul stellen, voor de waterbalans, voor $r=R_S$ levert nu

$$ Q_0 = \pi R_S^2 N $$

Of

$$ R_S = \sqrt{\frac{Q_0}{\pi N}}$$

Dit volgde natuurlijk al direct uit de waterbalans. Dit is ookk de reden dat de ligging van de waterscheiding voor variabele en constante pakketdikte dezelfde is.

BLom axiaal-symmetrisch

De situatie die Verruijt voor ogen heeft is dezelfde als die van Blom met uitzondering van wat er gebeurt buiten radius $R$. Bij Verruijt is radius $R$ een vaste rand met verlaging nul en waarbinnen de situatie kan worden opgevat als het gebied met drooggevallen sloten waarbinnen het neerslagoverschot niet meer wordt gedraineerd, maar richting de put stroomt. Bij Blom markeert de radius $R$ ook het gebied zonder sloten, maar stroomt er ook buiten deze radius nog water in de richting van de put. De grondwaterstroom $Q_R$ is zowel bij Verruijt als bij Blom gelijk aan $Q_0 - \pi R^2 N$. De grondwaterstroming over deze rand is bij Verruijt afkomstig van de vaste rand op $r=R$ en bij Blom komt die van grotere afstanden uit verminderde slootafvoer. Er is bij Blom op $r>R$ dus nog steeds verlaging in tegenstelling to bij Verruijt. Om dezelfde reden als hiervoor is de verlaging op $r=R$ gelijk aan $s_R = N c $.

Verlaging voor $r \le R$

De verlaging voor $r \le R$ is dezelfde als die bij Verruijt met dien verstande dat de velaging op $r=R$ niet nul is maar $s_R = N c $

Verlaging voor $r \ge R$ (De Glee (1930))

Bij Blom vallen de sloten droog binnen een straal gelijk aan $R$. De stroom $Q_R$ op $r=R$ volgt uit de waterbalans en is gelijk aan $Q_r = Q_0 - \pi R^2 N$. Het gebied voor $r>R$ fungeert als een semi-gespannen aquifer, die gevoed wordt uit verminderde drainage door de sloten, die gelijk is aan de verlaging ter plekke gedeeld door de drainageweerstand. De wiskundige oplossing voor de grondwaterstromign in deze situatie is volgens De Glee (1930); zij bestaat alleen voor constante pakketdikte:

$$ s_{r\ge R} = \frac{Q_R}{2 \pi kD}\cdot \frac{K_0 \left(\frac r \lambda\right)}{\frac R \lambda K_1\left(\frac R \lambda\right)} $$$$ Q(r) = Q_R \cdot \frac r R \cdot \frac{K_1\left(\frac r \lambda\right)}{K_1\left(\frac R \lambda\right)}$$

$Q_R$ volgt direct uit de waterbalans over het gebied gegeven door $r\le R$

$$ Q_R = Q_0 - \pi R^2 N$$

De verlaging op $r=R$ is gelijk aan $s_R$

$$ s_R = \frac{Q_R}{2 \pi kD}\cdot \frac{K_0 \left(\frac R \lambda\right)}{\frac R \lambda K_1\left(\frac R \lambda\right)} $$

De verlaging op $r=R$ is ook gelijk aan $s_R = N c$ omdat daar de verminderde drainage juist de resterende drainage nul maakt. Hiermee hebben we een vergelijking waarmee de radius kan worden berekend waarbinnen de sloten droog vallen

$$ N c = \frac{Q_R}{2 \pi kD}\cdot \frac{K_0 \left(\frac R \lambda\right)}{\frac R \lambda K_1\left(\frac R \lambda\right)} $$

We kunnen $R$ alleen iteratief berekenen, bijvoobeeld met de Newton methode. Deze methode maakt van bovenstaande expressie een functie van $y(R)$ waarvan $R$ het nulpunt is dat moet worden gevonden.

$$ y(R) = - N c + \frac{Q_R}{2 \pi kD}\cdot \frac{K_0 \left(\frac R \lambda\right)}{\frac R \lambda K_1\left(\frac R \lambda\right)} $$

De Newton methode vindt het nulpunt iteratief als volgt

$$R_{n+1} = R_N - \frac{y(R)}{y \prime{R}}$$

waarin $y \prime(R)$ de afgeleide is van $y(R)$ naar $R$ en $n$ de betreffende iteratie.

Deze methode is eenvoudig, maar vergt wel de afgeleide van $y(R)$. Deze kan zowel numeriek met $y(R)$ worden berekend als anlytisch. Voor dit laatste moeten we $y(R)$ differentiëren naar $R$.

Voor de afgeleide van $y$ naar $R$ hebben we de afgeleide van de factor met de besselfuncties nodig

Gebruik makend van $\frac{d K_0 z}{dz} = -K_1 z$ en $\frac{d (z K_1 z)}{dz} = -z K_0 z$ volgt na enige uitwerking

$$\frac{d}{dr}\left(\frac{K_0\left(\frac R \lambda\right)}{\frac R \lambda K_1\left(\frac R \lambda\right)}\right) = \frac 1 R \left( \frac{K_0^2 \left(\frac R \lambda\right)}{K_1^2\left(\frac R \lambda\right)} - 1\right) $$

waarmee uiteindelijk na invoegen in de gehele uitdrukking voor $dy(R)/dr$

$$\frac{dy(R)}{dr} = -\frac{\lambda N}{kD} \cdot \frac{K_0\left(\frac R \lambda\right)}{K_1\left(\frac R \lambda\right)} - \left(\frac{Q_0 / R}{2 \pi kD} - \frac{N R} {2 kD} \right) \cdot\left( 1 -\frac{K_0^2\left(\frac R \lambda\right)}{K_1^2\left(\frac R \lambda\right)} \right)$$

Deze functie $y(R)$ en $dy(R)/dR=y'(R)$ worden hierna geïmplementeerd en meegenomen in de Newton methoden om $R$ iteratief te bepalen. We zullen daarbij de analytische bepaalde afgeleide controleren met de numeriek bepaalde afgeleide.

Wanneer we $R$ kennen, kunnen we zowel de verlaging voor $r\le R$ als die voor $r\ge R$ berekenen.

Implementatie Verruijt en Blom eendimensionaal

Voor 1D berekening van de verlaging of de stijghoogte bij Verruijt moeten we afstand $L$ opgeven waar de verlaging nu is. Voor BLom is dit de afstand waarop de verlaging gelijk is aan $s_L = N c$. Bij Blom kunnen we die bij gegeven onttrekking en neerslagoverschot berekenen mits $kD$ en drainageweerstand $c$ bekend zijn.

Voorbeeld Verruijt 1D, vaste Q0, variërende L

Voorbeeld Verruijt 1D; vaste L, variërende Q0

De onttrekking en, daarmee de verlaging is stapsgewijs vergroot. Hierdoor komt de waterscheiding steeds verder weg te liggen tot deze de vaste rand $L$ bereikt.

Bij verruijt is buiten de vaste rand geen stroming en dus ook geen verlaging.

Bij grote onttrekking is de daling van de grondwaterstand in het eeste plaatje is een stuk groter dan de verlaging in het tweede plaatje. Dit is het gevolg van de afname van de dikte van het watervoerende pakket waar in het eerste plaatje wel en in het tweede plaatje geen rekening mee is gehouden.

Voorbeeld Blom 1D, variërende Q0

Bij Blom wordt de $L$ berekend, zodanig dat de verlaging op $x=L$ gelijk is aan $N c$.

In dit voorbeeld worden de grondwatertstand en de verlaging berekend voor verschillende waarden van $Q0$. Voor elke Q0 wordt de afstand $L$ berekend waarbinnen de sloten droogvallen. Op deze afstand is de verlaging gelijk aan $N c$, waardoor op afstand $L$ de sloot juist droogvalt (of beter: de sloot daar net niet meer draineert).

We zien verder dat de verlaging in het eerste plaatje, met variabele pakketdikte voor gotere onttrekkingen groter is dan die in het tweede plaatje, voor de situatie met vaste pakketdikte.

Voor de situatie met variabele pakketdikte zou de stijghoogte op $x > L x$ nog iets gecorrigeerd kunnen worden voor de in werkelijkheid afnemende dikte. Dit effect is echter zo klein dat het verschil in de aansluiting op het aangegeven punt, dus $x=L$ in de grafiek niet is te zien. Deze correctie kan eenvoudig worden verwaarloosd in de praktijk.In het onderhavige geval is dit een correctie van de $\lambda van $h/H = (D - s_L) / D \approx 19.5 / 20 \approx 0.98$ op de gebruikte waarde van $\lambda$, dus verwaarloosbaar.

Verruijt en BLom, Axiaal symmetrisch

Implementatie

Voorbeelden Verruijt, axiaal symmetrisch

Voorbeeld Blom, axiaal symmetrisch

Showing the progress of the Newton iterations to find the point R in y(R)=0

De afstand van de put waarop de verlaging precies gelijk is aan Nc, het critierium voor juist droogvallen van de sloten wordt iteratief berekend met de methode van Newton. Het voorschrijden van dit iteatieproces wordt hieronder grafisch weergegeven.

Voor de iteraties is de afgeleide nodig van de fuctie $y(R)$ zie boven. De tweede grafiek toont de afgeleide, zowel analytisch als numeriek berekent ter controle.

Ruimtelijk modelleren van verlagingen met de drainageweerstand (experimenteel, het is de vraag of dit werkt)

Een oppervlaktewatrstelsel in een gebied kan op regionale schaal mathematisch of modeltechnisch in de berekeningen of in de modellering worden verdiscondeerd door deze ruimtelijk om te zetten in een drainageweerstand. De drainageweerstand is, per definitie, de gemiddelde grondwaterstand minus het oppervlaktewaterpeil, gedeeld door het gemiddelde neerslagoverschot. Immers $\frac{h_{gw}-h_{ow}}{c}=N$, dus [m/d]=[m/d]. In feite is deze defitie van de drainageweerstand een regionale meting ervan, en hij kan dus worden opgevat als een eigenschap van het concrete, zelfs fysische systeem van grond- en oppervlaktewater op sub-regionale schaal, waarin het oppervlaktewater met zijn individuele sloten tot de regionale drainageweerstand zijn opgeschaald en dus uit de vergelijking zijn geëlimineerd.

De drainageweerstand kan ruimtelijk worden vastgesteld, bijvoorbeeld door een grondwatermodel door te rekenen zonder neerslagoverschot en daarna nog een keer, maar dan met het gemiddelde neerslagoverschot, en dit verschil te delen door het neerslagoverschot. Dit is dus een exercitie die direct met elk grondwatermodel kan worden uitgevoerd.

De volgende stap is dan om in hetzelfde of in een ander model al het oppervlaktewater te vervangen door deze drainageweerstand. Met dit model kan dan heel eenvoudig het effect worden berekend van ingrepen zoals onttrekkingen, bemalingen en aanpassingen van oppervlaktwaterpeilen. Het effect van zulke aanpassingen van het oppervlaktewaterpeil kan worden berekend door uitsluitend de in peil te vanderen oppervlaktewateren in te voeren, samen met hun peilverandering.

Voor deze wijze van modellering zijn de absolute peilen niet van belang. Zelfs voor de berekening van de drainageweerstand is het niet nodig om de absolute hoogte van het oppervlaktewater te kennen. Immers het effect van het veranderen van de grootte van het neerslagoverschot op de grondwatersanden is een zuivere superpositie, waarvan het resultaat geheel onafhankelijk is van de (stationaire) beginsituatie. Het enige gege ens die we voor een dergelijke analyse nodig hebben zijn de lagen van het model met zijn bodemconstanten en, aan de top, de ruimtelijk variërende drainageweeerstand.

Om de drainageweerstand te kunnen bepalen is het zaaks de gemiddelde grondwaterstand ten opzichte van het lokale opervlaktewater te kennen. Het peil in het oppervlaktewater, lokaal of niet, kan daarbij op nul worden gezet. Het berekende ruimtelijke grondwaterstandsverloop varieert, zodat een slag moet worden gemaakt naar het gemiddelde verloop dat een goede maat oplevert voor de drainageweeerstand. Een mogelijkheid is om afzonderijke gebieden in het model met oppervlakgewater te selecteren waarin de grondwaterstand een opbolling vertoont als gevolg van de opgelegde voeding. Met GIS kan voor ek van deze deelgebieden de gemiddelde grondwaterstand worden berekend. Dit is dan de maat voor de drainageweerstnad binnen elk van de afzonderlijke gebieden met opbolling.

Het vinden van een gebied met opbolling is een GIS-exercitie. Voor de uitvoering zijn nog wel een paar hobbels te nemen. Als het regionale gebied is op te delen in gebieden met zuivere opbolling, dat wil zeggen, dat de tweede afgeleide van de grodnwaterstand in alle richtingen positief is, is er weinig tegen deze analyze in te brengen. Dat is echter anders in overgangsgebieden, dit zijn gebieden die van een hoog gebied zonder oppervlaktewater overgaan in een laag gebied met veel oppervlaktewater. In het overgangsgebied prikken beken vanuit het lage gebied stroomopwaarts het hoge gebied in. In deze overgang zijn de gebiedsdelen met meer beken open in de richting van het hoge gebied en gaat de kromming van de stijghoogte van het lagere gebied naadloos over in die van het hogere gebied. De vraag is dus hoe we in deze open overgangsgebieden op een zonvolle manier een drainageweerstand kunnen definiëren, oftewel een zinvol gemiddelde van de opbolling kunnen bepalen, zo mogeijk door middel van een standaard GIS-bewerking. Vermoedelijk ligt de sleutel in het vinden van zadelpunten in het grondwateroppervlak, want hierin gaat de tweede afgeleide van het grondwateroppervlak in bepaalde richtingen door nul. Een eigenschap van een opbollingsvlak is dat de tweede afgeleide van het opbollingsvlak groter of gelijk aan nul is in alle richtingen. Waar dat niet het geval is, hebben we te maken met een overgang in de drainageweeerstand.

Het verhaal in de vorige paragraaf is voor de praktijk veel te ingewikkeld. Maar omdat de tweede afgeleide van het grondwatervlak nooit voldoende nauwkeurig kan worden bepaald, kan zo'n aanpak ook niet leiden tot een zinvol en stabiel resultaat. Vandaar dat hierna een andere benadering wordt voorgesteld. Die echter nog niet is uitgetest en dus zuiver experimenteel is.

Voorstel werken met drainageweerstand (experimenteel)

De voeding van een gebied kan afgeleid worden uit de bolling (tweede afgeleide naar $x$ en $y$) van het grondwatervlak, immers $$\frac N {kD} = \frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial y^2} =\nabla^2 h$$

Voldoet de grondwaterstand aan deze partiële differentiaalvergelijking, dan is de voeding overal gelijk, met de waarde $N$. Als we nu een model zoeken waarbij de lek bij opgelegde werkelijke grodnwaterstand overal gelijk aan $N$ is. Alleen eisen we dit nu niet door voeding van het watervoerende pakket vanuit de lucht, maar als gevolg van de opgelegde gemiddelde grondwaterstanden die via een weerstandbiedende deklaag met (drainage)weerstand $c_{dr}$, in verbinding staat met het slootpeil"

$$ \frac {h - h_0} {c_{dr}} = N $$

Met $h_0$ de grondwaterstand zonder voeding. Vanwege de geldigheid van het superpositiebeginsel, mag $h_0$ gelijk aan nul worden genomen.

De ruimtelijk variabele drainageweerstnad is per definitie gelijk aan

$$ c_{dr}(x, y) = \frac {h - h_0} N$$

De drainageweerstand kan zo worden verkregen als zijnde $h - h_0$ uit ons testmodel bij constante voeding gelijk aan $N = 1$ m/d.

Met deze weerstand kunnen we vervolgens de effecten van ingrepen berekenen, waarmee dan automatisch de invloed van het oppervlaktewater is meegenomen.

Voorbeeld

In het navolgende voorbeeld beschouwen we een gebied met een gemiddelde, ruimtelijk uniform neerslagoverschot en oppervlaktewater. De grondwaterstand is overal met en model uitgerekend.

Hieruit bepalen we de drainageweerstand, $c_{dr} = \frac {h - h_0} N$$

Om te laten zien dat het niet uitmaakt wat de waterstanden in het oppervlaktewater zijn, doen we hetzelfde met alle waterstanden op nul gezet.

Dan draaien we een ander model, dat hetzelfde is als het voorgaande, maar dat in plaats van oppervlatewater bevat een weerstandbiedende laag aan de top bevat met de drainageweerstand. De stijghoogte (vervangend oppervlaktewaterpeil) boven de afdekkende laag is overal nul, het neerslagoverschot is ook nul, maar met de hierboven de berekende gemiddelde grondwatestand overal opgelegd, is de stroming nu gelijk aan die in het eerste model en door de draiageweerstand is de voeding nu overal gelijk aan N, en dus is de berekende grondwaterstand nu ook overal gelijk aan de in het vorige model beekende grondwaterstand. In dit alternatieve model kunnen vervolgens maatregelen worden opgelegd en het effect daarvan berekend. Dat kan gewoon met all oppervlaktewaterpeilen gelijk aan nul. Het resultaat is een verlagingenbeeld dat exact overeenstemt met de drainageweerstand en goed met de verlagingen die we zouden verkrijgen met het model waarin het oppervlakterwater exact was gemodelleerd.

Deze test moet worden uitgevoerd met een numeriek model zoals Modflow, aangezien TTIM geen ruimtelijk variabele drainageweerstand kan moddelleren.

Het bijzondere is dat voor deze exercitie geen GIS nodig is. Bovendien kan deze alternatieve modellering met drainageweerstand met een zeer eenvoudig model worden uitgevoerd, zonder oppervlaktewater; immers het effect daarvan is verwerkt in de drainageweerstand die ruimtelijk varieert.

Dit is nog uit te proberen.

Hoe lang duurt het voordat de verlagingen als stationair kan worden beschouwd?

Bij onttrekking met een put in een oneindig uitgestrekt watervoerend paket ontstaat nooit een stationaire eindsituatie omdat er geen terugkoppeling is tussen verlagen en infiltratie.

Verlagingen kunnen alleen stationair worden wanneer ergens in het watervoernde pakket randen zijn met vaste stijghoogte, danwel er een vaste stijghoogte is in een bovenliggend pakket waarmee hydraulische verbinding bestaat via lek door een semi-doorlatende laag. De essentie is dat door de verlaging een stroming wordt opgewekt die op termijn de onttrekking volledig kan compenseren.

De belangrijkste situatie waarbij stationaire stoming optreedt is meestal die met oppervlaktewater verderop in het gebied of een semi-permanente toplaag waarboven het oppervlaktewaterpeil ruimtelijk wordt beheerd.

Voor een put die onttrekt uit een pakket met semi-gespannen water geldt de formule van Hantush (1955)

$$ s(r, t) = \frac{Q}{4 \pi kD} W_h \left(u, \rho\right), \,\,\,\,with\,\,\,\, u = \frac{r^2 S}{4 kD t},\,\,\,\,and\,\,\,\,\rho=\frac{r}{\lambda}$$

Wanneer deze verlaging wordt uitgezet op logartmische tijdschaal blijkt er een buigpunt te zijn halverwege de eindverlaging, waar de verlaging steeds trager naar toe loopt. Het tijdstip van dit buigpunt is een geschikt moment om op praktische wijze de tijd te bepalen waarna de verlaging als stationair mag worden beschouwd.

Het buigpunt ligt op

$$ u = \frac \rho 2 = \frac{r}{2 \lambda}$$

Om het snijpunt van de raaklijn door het buigpunt en de horizontale lijn van de eindverlaging te bepalen, is het mathematisch handiger om direct met de de parameter op de horizontale as te werken

$$ \eta = \ln \left(\frac 1 u \right) = -\ln u $$

De tangent of the raaklijn door de het buigpunt op de as $1/u$ is

$$ \left(\frac{\partial W_h}{\partial \eta} \right) = e^{-\rho}$$

De stationaire eindwaarde is

$$W_h(0, \rho) = 2 K_0(\rho)$$

Waar de raaklijn door het buigpunt de lijn met de eindverlaging snijdt is de verlaging halverwege de eindverlaging behorend bij de stationaire situatie. Het punt waarop de verlaging stationair is geworden kan praktisch worden gedefinieerd als het punt waarop de raaklijn door het buigpunt de horizontale lijn met de eindverlaging snijdt.

$$ \eta = e^{\rho} K_0(\rho) + \eta_{Bp},\,\,\,\,with\,\,\,\,\eta = -\ln(u)$$$$ \eta_{Bp} = -\ln \left(u_{Bp} \right),\,\,\,\, with\,\,\,\,u_Bp = \frac{\rho}{2}$$$$\eta - \eta_Bp = e^{\rho} K_0(\rho)$$$$\frac{e^\eta}{e^{\eta_{Bp}}} = e^{e^{\rho} K_0(\rho)}$$$$\frac{u_{Bp}}{u} = e^{e^{\rho} K_0(\rho)}$$

Het stijdtip (de $u$-waarde) waarop de verlaging als stationair mag worden beschouwd volgt dus direct uit de $u$-waarde van het buigpunt $u_{Bp}$ en $\rho=r/(2 \lambda)$:

$$ u = u_{Bp} e^{-e^{\rho} K_0(\rho)}$$

In het voorbeeld hieronder wordt de Hantush functie getoond versus $1/u$ op logaritmische schaal. Het eerste deel van de dubbele buiging start wanneer de invloedsgrens van de onttrekking het waarnemingsput bereikt en het tweede deel wanneer de verlaging wordt afgeremd door de opgewekte lek. De eindverlaging is gelijk aan de stationaire verlaging, zodat, met $t\rightarrow\infty,\,\,u\rightarrow0$:

$$ W_h(0, \rho) = 2 K_0(\rho)$$

Het buigpunt light op halve hoogte, dus bij $W_h(u_{bp}) = K_0(\rho)$. Voor het buigpunt geldt bovendien

$$u_{bp} = \frac{\rho}{2} = \frac{r}{2 \lambda} $$

De berekende buigpunten worden met een dikke stip in de grafiek aangegeven.

Voorts is in de grafieken aangegeven waar de raaklijn door het buigpunt en waar de eindverlaging snijdt (stippelijn met snijpunt).

Het snijpunt van deze raaklijn met de eindverlaging ligt op

$$ u = u_{Bp} e^{-e^{\rho} K_0(\rho)}$$

Deze gemakkelijk te berekenen waarde van $u$ is een praktische waarde voor het definieren van het moment waarop de eindverlaging is bereikt kan worden beschouwd.

$$ u = \frac{\rho}{2} e^{-{e^\rho} K_0(\rho)}$$

Voor een reële situatie (met het argument in de Theis-formule (geen lek) links staat en de invloed van de lek rechts):

$$u = \frac{r^2 S}{4 kD t}=\frac{r}{2 \lambda} e^{-{e^{\frac{r}{\lambda}} K_0\left(\frac{r}{\lambda}\right)}}$$

of

$$ u = \frac{r^2 S}{4 kD t}=\frac 1 2 \rho e^{-{e^{\rho} K_0(\rho)}}$$

Het voorbeeld toont het verloop van de verlaging in de tijd, met de tijd op logarithmische schaal voor verschillende waarden voor $\lambda$. Op het buigpunt is de verlaging steeds halverwege de stationaire eindwaarde. Een redelijk critirum voor het bereiken van de stationaire eindwaarde lijkt het snijpunt van de raaklijn door het buigpunt met de horizontale lijn die de eindverlaging weergeeft. Dit punt is steeds met een grafische punt weergegeven in de grafiek.

Onvolkomen putten

Een onvolkomen put veroorzaakt op kortere afstanden van de put een afwijking van de verlaging ten opzichte van de situatie met een put met zijn filter over de gehele dikte van de watervoerende laag. Onvolkomen filters zijn echter de gangbare situatie, met name in dikkere pakketten.

We kunnen de invloed van de onvokomenheid van het putfilter verwerken met een correctie van Hantush. Deze correctie kan gewoon worden meegenomen in de analytische oplossing.

Overigens is het effect van de onvolkomenheid van het filter uitgewerkt op afstanden groter dan circa anderhalve pakketdikte en daarmee in de regel niet van belang voor de vergunningverlening.

Nog uitwerken (kan eenvoudig met formule hiervoor van Hantush (1956?) of eventueel van Huisman (1970)).

Vertraagde nalevering (Delayed yield)

Pompproeven die als basis voor aanvragen van vergunningen worden uitgevoerd, zijn vaak zo kort dat alleen de verlaging in het bepompte pakket kan worden bepaald, en geen informatie wordt verkregen overhet na langere tijd nazakken van het freatisch vlak boven de scheidende laag. Dit freatische vlak zakt later na omdat het pompen lek opwekt door de scheidende laag op het watervoerende pakket. Dit nazakken is vertraagd ten opzichte van de verlaging die in het watervoerend pakket optreedt, wat het gevolg is van de bergingscoëfficiënt die voor het freatisch pakket veel hoger is dan die van het semi-gespannen watervoerende pakket. De voortoets zou altijd moeten eisen dat dit effect van nazakken wordt meegenomen in de analyse omdat de eindwaarde die ogenschijnlijk na enkele dagen al is bereikt in het watervoerende pakket niet de eindwaarde van de verlaging in het freatisch pakket is en derhalve misleidend en vaak ten onrechte wordt gebruikt als eindverlaging.

$$ s = \frac Q {4 \pi kD} W(u) = \frac Q {2 \pi kD} K_0 \left(\frac r \lambda \right) $$$$ W(u) = 2 K_0 \left(\frac r \lambda \right) $$$$ \ln \left( \frac {2.25 kD t}{r^2 S_y} \right) = 2 K_0 \left(\frac r \lambda \right) $$$$t = \frac{r^s S_y} {2.25 kD} e^{2 K_0 \left( \frac r \lambda \right)}$$

De algemene situatie met vertraagde navlevering kan voor het semi-gespannen pakket worden benaderd door de aanvankeijke verlaging volgens Hantush met de elastische bergingscoëfficiënt S te laten vervolgen met de verlaging volgens Theis met de freatische bergingscoëfficiënt. Het aansluitpunt van de twee curves is het punt waarop de freatische verlaging voglens Theis die van Hantush inhaalt.Vanaf dat moment zijn de verlaging in het ondiepe en diepe pakket geijk, en lopen samen verder op. In de werkelijkheid ontstaat een vloeiende overgang van de Hantush verlaging naar die volgens Theis, die met een tweelagen model kan worden berekend. Dat is wat ingewikkelder en niet noodzakelijk voor een adequaat beeld van wat er in de werkelijkheid gebeurt. Voor de praktijk kan de hier gevolgde methode als benadering worden gebruikt. Hij laat in elk geval zien hoe lang men een pompproef doen om zeker te weten dat er in later stadium geen substantiële daling van het freatisch vlak optreedt.

In het voorbeeld hieronder met de verlaging op drie afstanden in het watervoerend pakket zijn de getrokken lijnen berekend met de formule van Hantush (elastische bergingscoëfficiënt) voor semi-gespannen water. De stippellijnen zijn berekend met de formule van Theis (freatische bergingscoëfficiënt). Na verloop van tijd halen de lijnen volgen Theis die volgens Hantush in. Vanaf dat moment gaan de verlaging van het freatische water en die van het semi-gespannen water gelijk op. De oplopende lijnen vanaf het intersectiepunt is de nazakking, die zowel in het freatische als in het sem-gespannen pakket plaats vindt.

Worst case = "Je mag het niet missen, wanneer de voorgenomen ingreep serieuze gevolgen heeft"

Hoe vullen we dat in? Nog uit te werken. Mogelijk niet eenduidig.

  1. Rekenen tot wanneer de verlaging statoinair is geworden.
  2. De nazakking niet vergeten.
  3. Het "intrekgebied" volgens de formule van Verruijt is niet het gebied waarin het grondwater wordt verlaagd. Die wordt bepaald tot de afstand waarop de verlaging nul is, een vaste rand. De verlaging wordt gegeven door Dupuit en niet door Verruijt, tenzij de voeding bij Verruijt het gevolg is van afgenomen drainage door droogvallende sloten.
  4. De juiste aanpak in een gebied met gebiedswijze drainage door sloten is die volgens Blom. Hierbij worden de sloten in een gebied rond de put geheel droog getrokken en komt de drainage door die sloten te vervallen en ten goede aan de onttrekking. De rand van dit gebied is die waarbij de verlaging gelijk is aan de drainageweerstand maal het gedraineerde neerslagoverschot. Daarbuiten fungeren de nog niet geheel drooggevallen sloten als voeding die zich gedraagt als voeding vanuit een bovengelegen laag met vast peil en een drainageweerstand.
  5. De formules van Verruijt en Blom zijn zowel ééndimensionaal als axiaal-symmetrisch toe te passen.

Grondwaterbalans. Hoe hiermee omgaan in de passende beoordeling?

Dit gaat onder andere om indirect, via het grondwater, onttrekken aan beken met beperkte afvoer, die in de zomer kwetsbaar zijn. Dit speelt in het algemeen op hogere gronden waar geen externe aanvoer naar het oppervlaktewater kan plaats vinden, dus waarin al het oppervlaktewater afkomstig is van toestromend grondwater.

De onttrekking eigent zich water toe dat anders een andere bestemming zou vinden, zoals voeding van beken en of van de vegetatie. Er zijn derhavle atijd gevolgen van een onttrekking voor de waterbalans van het gebied waarop de ontrekking invloed uitoefent. Men kan berekenen hoeveel indirect via het grondwater aan een naburige beek wordt onttrokken, repectievelijk aan voeding wordt ontnomen. Met name kleinere beken zijn gevoelig voor aanpalende grondwateronttrekkingen. Regionaal kan het cumlatieve effect grote negatieve gevolgen hebben voor de beekafvoer en daarmee voor de ecologie en andere gebruiksdoelen. In de literatuur zijn vele voorbeelden te vinden waarin beken en hele rivieren droog zijn gevallen door onttrekking van grondwater.

GxG kaart. Welke rol kan die spelen bij de passende beoordeling?

Kan hiermee Vlaanderen worden verdeeld in grondwatergebiedstypen: Infiltratiegebied, intermediar en kwelgebied? Zulke kaareten geven aan wat er als natuurlijke variatie in de grondwaterstand beschouwd kan worden.

Er kan voor de beoordeling van vergunningaanvragen rekening worden gehouden met het type grondwatersysteem: infiltratiegebied (hoog gelegen met diepe grondwaterstand, weinig of geen drainagemiddelen of natuurlijke afvoer) en een kwelgebied (relatief laag gelegen ten opzichte van de omgeving, met hoge grondwaterstanden en veel drainagemiddelen), danwel een overgangsgebied tussen deze twee uitersten in. Naar verwachting zullen toekomstige GxG kaarten deze gebiedsindeling weerspiegelen. De gevolgen van een onttrekking in de verschillende zones zullen ook verschillen. Zo zal een onttrekking in een infiltratiegebied ogenschijnlijk nabij nauwelijks invloed hebben maar zal de invloed verderop manifesteren in het droogvallen van koppen van beken, en van sloten in het overgangsgebied. Deze kunnen daar heel gevoelig voor zijn, terwijl het een op een aantonen van deze effecten aan de hand van metingen wordt bemoeilijkt door de grote natuurlije variatie van de afvoer van de beekkoppen en de van nature variënde slootpeilen in de overganszone. De reikwijdte van de invloed van de onttrekking in een infiltratiegebied is in het algemeen groot doordat de randen die de invloed dempen ver weg liggen.

De situatie bij een onttrekking in een kwelgebied met veel oppervlaktewater, dat stevig gevoed wordt door kwel uit de verre omgeving en (dus uit een groot gebied) is omgekeerd. De reikwijdte is door het vele dempende oppervlaktewater in de naaste omgeving beperkt, zowel wat de verlaging betreft als wat betreft de invloed op het oppervlaktewater.

De situatie in de overgangszone ligt hier tussenin. Beken en sloten tegen het hogere gebied aan zijn zeer gevoelig voor de onttrekking, voor die in de richting het kwelgebied is dit juist minder het geval. De invloed van de onttrekking en zijn reikwijdte is daarom asymmetrisch, hij is groter richting het infiltratiegebied en kleiner richting het kwelgebied.

Deze nuances zijn analytisch moeilijk te benaderen. Een model heeft hier echter geen moeite mee. Analytisch zou men de invloedsradius kunnen kiezen afhankelijk van het gebiedstype, af te leide uit de GxG kaarten. Een methode hiervoor zou moeten worden uitgewerkt en getoetst aan een aantal representatieve situaties.

Invloedstraal (radius of influence)

Om de in de tijd toenemende reikwijdte van een onttrekking aan te geven wordt wel gesproken van invloedsgrens of "radius of influence". Deze "radius of influence" bij een continue onttrekking in een zich oneindig uitstrekkend watervoerend pakket zonder lek kan direct worden afgeleid uit de verlaging volgens Theis, die in dit geval van toepassing is. De verlaging volgens Theis kan voor voldoende grote tijd worden benaderd met de vereenvoudigde formule met daarin de logarithme in plaats van de well function:

$$ s = \frac Q {4 \pi kD} \ln\left(\frac{2.25 kD t}{r^2 S}\right) $$

Voor het nog niet beïnvloede gebied geldt dat de verlaging nog nul is. Dit is hier (bij benadering) het geval wanneer het argument onder de logaritme gelijk is aan 1:

$$ \frac{2.25 kD t}{r^2 S} = 1 $$

zodat daar een "radius of influence" uit volgt, die gelijk is aan

$$ r = \sqrt{\frac {2.25 kD t}{S}} $$

Deze formule is een eenvoudige doch uitermate praktisch uitdrukking om de radius te bepalen waarbuiten de invloed van de onttrekking nog niet bemerkbaar is. Het is een benadering, omdat de gebruikte formule van de verlaging een benadering is, maar het is een praktische benadering die op de aangegeven wijze objectief is afgeleid. Hieronder volgt een voorbeeld waarin de werking van deze radius inzichtelijk wordt gemaakt.

Exacte reikwijdte voor de Theis situatie

De exacte tijdsafhankelijke reikwijdte (voor zover die in de praktijk bepaalbaar is) bij een freatische verlaging van bijv. 5 cm, vergt een wat nauwkeuriger berekening dan de hiervoor afgeleide eenvoudige formule. Daarvoor hebben we wel de inverse nodig van $W(u)$, dus de $u$ die hoort bij bekende $W(u)$.

$$ s = \frac Q {4 \pi kD} W(u) $$

en dus $$ W(u) = 4 \pi kD \frac s Q $$

waarbij $s$ gegeven is net als $Q$. Om de reikwijdte te berekenen hebben we $u$ nodig

$$ u = W^{-1}(u) $$

Voor de inverse van W(..) (de exponential integral) bestaat geen kant en klare analytische uitdrukking. Ergo de inverse moet numeriek worden uitgerekend. Dit kan snel met Newton's methode of nog sneller met Halley's method (zie wikipedia). Echter het is aan te bevelen om in plaats van W(u) ln(W(u) als functie te nemen en ln(u) als argument. We beperken ons vanwege de afleiding dan wel tot Newton's method.

Oplosseen voor

$$ F(z, A) = ln(W(u)) - ln(A) = 0$$

met $z = ln(u)$, zodat, met $u = e^z$ volgt dat $\frac{d e^z}{dz} = e^z=u$

$$ F1(z) = \frac{dF(z, A)}{dz} $$

$$ \frac{dF(z, A)}{dz} = \frac{d(\ln(W(u)) - \ln(A))}{dz}) = \frac{1}{W(u)}\frac{dW(u)}{du}\frac{du}{dz}= \frac{1}{W(u)}\frac{-e^{-u}}{u}\frac{de^{z}}{dz}= \frac{1}{W(u)}\frac{-e^{-u}}{u} u=\frac{-e^{-u}}{W(u)}$$

$$ F2(z) = -\frac d {dz} \left(\frac{e^{-u}}{W(u)}\right) = W^{-2}(u)\frac{e^{-u}}{u}\frac{du}{dz} e^{-u} - W^{-1}(u)e^{-u}\frac{du}{dz} =\frac{u e^{-u}}{W(u)}-\frac{e^{-2u}}{W^2(u)}$$

Newtons's methode om het nuplpunt te vinden $$ z_{n+1} = z_n - \frac{F(z_n, A)}{F1(z_n)} $$

en, uiteindelijk de terugvertaling naar $u$ $$ u = \ln(z) $$

Halley's method convergeert nog sneller dan die van Newton en is gedefineerd als (Wikipedia)

$$ z_{n+1} = z_n - \frac{F(z_n, A)}{F1(z_n)} \left[1 - \frac{F(z_n, A)}{F1(z_n)} \times \frac{F2(z_n)}{2 F1(z_n)}\right]^{-1} $$

Halley's methode is als optie ingebakken in de implementatie van de funktie beneden. Voor Halley's methode moet $A = W(u)$ < 20, wat in de praktijk altijd zo is.

Verificatie van F(z, wu) en de afgeleiden F1(z) en F2(z)

De verificatie laat het volgende grafisch zien. De blauwe streepjeslijn is de well-function van Theis versus $u$ (in plaats van de gebruikerlijke $1/u$) op logschaal. De rode punt is het punt waarvan $W(u)$ is gegeven en waarvan $u$ wordt gezocht. De getrokken oranje lijn is dezefde grafiek verminderd met de waarde W(u), het is de functie $y(u)$. De rode punt wordt nu de blauwe punt waarvoor $y(u)=0$. De Newton methode zoekt nu iteratief waar de functie $y(u)$ door nul gaat. Tevens zijn de eerste en de tweede afgeleide van de functie getekend, zowel analytische als numeriek berekend, ter controle.

De curves laten zien dat de afleidingen juist zijn en kunnen worden toegepast in de Newton methode en voor nog wat snellere convergentie in de methode van Halley.

Implentatie van de inverse van de well-function van Theis (exp1(u) <-> exp1_inv(W))

Met het voorgaande kan een functie worden gedefinieerd die de inverse is van de well function van Theis, die we hier exp_inv(..) noemen met als argument de waarde van W($u$) en als uitkomst die van $u$.

Tegelijk wordt de uitkomst gecontroleerd door $u$ weer in de functie exp1(..) in te vullen.

Toepassing vergelijking van de eenvoudige en de meer complexe berekening van de invloedsstraal

Hieronder wordt de invloedsstraal berekend met de eenvoudige methode $r=\sqrt{\frac{2.25 kD t}{S}}$ an met de exacte verlaging $s=0.05$ m, waarvoor de boven geïmplementeerde funtie exp1_inf(s=0.005) wordt gebruikt. Het verschil blijkt voor de praktijk verwaarloosbaar te zijn.

Berekening van de verlaging, ten gevolge van een willekeurig in de tijd verlopende onttrekking, met convolutie

Convolutie is een algemene method waarmee het dynamische verloop van dde grondwaterstand of verlaging kan worden berekend die volgt uit een in de tijd variërende onttrekking of andere stress zoals neerslagoverschot. Convolutie is in feite een slimme manier van superpositie, waarmee in eenslag een hele tijdreeks kan worden berekend.

De stanaard oplossingen voor het niet-stationaire verloop van de verlaging door een stationaire onttrekking zijn op te vatten als stapresponses wanneer de onttrekking gelijk aan 1 wordt genomen. Bijvoorbeeld de verlaging volgens Theis. De staprespons volgens theis is dan

$$SR = \frac{1}{4 \pi kD} W_t(u) = \frac 1 {4 \pi kD}\intop_u^\infty \frac{e^{-y}}{y} dy,\,\,\,u=\frac{r^2 S}{4 kD t}$$

Hieruit volgt de zogenoemde impulsrespons door differentiatie

$$ IR = \frac{\partial}{\partial t} SR = \frac Q {4 \pi kD} \frac{e^{-u}}{t}$$

De impulsrespons heeft een groot theoretisch belang, maar is voor ons niet nodig. Belangrijker in de praktijk is de staprespons die het gevolg is van een continue onttrekking van 1 m3/d vanaf t=0 en de daaruit afgeleide blokrespons die de reactie is van de grondwaterstand op een plotseling onttrekking van 1 m3/d gedurende exact 1 tijdstap (meestal een dag)

$$ BR = \frac Q {4 \pi kD} \left[ W_t\left({u_t}\right) - W_t\left({u_{t-dt}}\right)\right] $$

De berekening van de verlaging door een willekeurig verlopend debiet geschiedt het gemakkelijkst met de functie lfilter(...) uit de modulescipy.signal.

$$s(r,t) = \frac{1}{4 \pi kD} \sum_{i=0}^\infty Q(t - \tau_i) BR(\tau_i) =lfilter \left(BR(\tau, 1, Q(t-\tau)) \right)$$

Het werken met de blokrespons is doorgaans het meest effectief. Het is hetzelfde als convolutie.

In het voorbeeld hieronder wordt de verlaging volgens Theis gegeven als formule (lijn) en zoals berekend met convolutie met de blokrespons (puntjes) Het blijkt dat de berekening met convolutie, gebruik makend van de blokrespons hetlzelfde resultaat geeft als de analytische oplossing van Theis. Dit moet natuurlijk ook zo zijn. We kunnen hiermee dus ook de verlaging volgens Theis berekenen bij een willekeurig in de tijd verlopende onttrekking.

Convolutie met de impulsrespons in plaats van de blokrespons

We kunnen de convolutie ook uitvoeren met de impulsrespons. Dit is de afgeleide van de staprespons en is de reactie van het grondwatersysteem op een pulsvorminge onttrekking. Dat wil zeggen een onttrekking der lengte $dt\rightarrow0$ ter grootte van een eenheid, dus een onttrekking van 1 m3 gedurende een tijdstap $dt$ die infinitesimaal klein is.

$$ s(t) = \intop_{\tau=0}^\infty Q(t - \tau)\frac{\partial SR}{\partial \tau} d\tau \approx\sum_{i=0}^\infty Q(t-i \Delta \tau)\frac{\partial SR}{\partial \tau}\Delta \tau$$

Waarbij de rechter uittdrukking de numerieke convolutie voorstelt met discrete tijdstappen van $\Delta \tau$, veelal van 1 dag.

Het probleem hiermee is dat de impuls respons zo snel verandert binnen het tijdsbestek van 1 dag, dat de berekening onnauwkeurig wordt en de uitkomst ver kan afwijken van wat hij moet zijn. We kunnen de uitkomst wel goed krijgen, maar dan moeten we de tijdstap veel korter kiezen dan 1 dag. Dit effect wordt hieronder gedemonsteerd, waar we de verlaging volgens Theis opnieuw berekekene, maar nu met zowel de blokrespons als de impulsrepons. Voor grote tijdstappen blijkt de impulsrespons veel te kleine waarden op te leveren; pas bij tijdstappen kleinder dan 0.03d komen beide berekeningen goed met elkaar overeen.

In de praktijk gebruiken we altijd de blokrespons, want die is exact en niet afhankelijk van de grootte van de tijdstap.

We kunnen nu, door convolutie met de blokrespons, nauwkeurig de verlaging bij elk willekeurig verloop van de onttrekking uittrekenen.

Enkele voorbeelden van convolutie met niet constante onttrekking

  1. Willekeurig verloop van de onttrekking
  2. Onttrekking neemt lineair af in de tijd, waarbij de onttrekking gehalveerd wordt over een tijdvak T
  3. Onttrekking neemt exponentiaal af in de tijd, waarbij de onttrekking halveeft over een tijdvak T

Welke onttrekking is nodig om de gewenste verlaging op een gegeven afstand te handhaven?

Bemalingen hebben als doel een verlaging overal binnen de bouwput te bereiken om zo droog te kunnen werken en bouwen. Er zal normaliter een voldoend hoog debiet worden ingezet om binnen een redelijke tijd de vereiste verlaging te behalen. Daarna zal het debiet geleidelijk worden verminderd omdat de verlaging bj gelijk debiet anders steeds verder zou oplopen. Deze geleidelijke vermindering van het debiet is dus noodzakelijk om niet nodeloos veel te pompen of teveel te verlagen.Kunnen we het verloop van het debiet gedurende de gehele duur van de bemaling berekenen?

Een bouwput bemaling kan schematisch worden opgevat als een uniforme peilverlaging die wordt verkregen door pompen binnen een een cirkel met straal $\rho$, de geschematiseerde rand van de bouwput. We kunnen ons het bemalen voorstellen als een schermvormige onttrekking met een totale omvang $Q_0$ m3/d, die gelijkmatig is verdeeld langs de cirkelvormige rand van de bouwput. Zoals eerder in dit document is aagetoond, is de verlaging binnen de bouwput dan vlak en bovendien gelijk aan de verlaging ter plaatse van de rand die zou worden veroorzaakt door bemaling met een enkele punt met hetzelfde totale debiet in het centrum van de bouwput. Voor de verdere analyse kunnen we dit als uitgangspunt nemen.

De analytische onttrekkingsformules die ons ter beschikking staan vergen alle een gegeven vast onttrekkingsdebiet, waarvan de verlaging dan het gevolg is. Uit de literatuur zijn wel analytische formules bekend die de verlaging opleggen en daarbij het in de tijd verlopende bemalingsdebiet uitrekenen, maar die zijn veel complexer dan die met opgelegd debiet. De betreffende formules zijn gegeven in de vorm van een integraal met daarin verschillende oscillerende Besselfuncties die moeilijk goed numeriek zijn te integreren.

Hieronder wordt een zeer eenvoudige procedure voorgesteld en uitgwerkt om de gewenste verlaging op geggeven afstand, de rand van de bouwput met radius $\rho$, in de tijd te handhaven. De methode is niet exact, maar wel een goede en praktische benadering.

De verlaging in de tijd in de situatie van Theis is

$$ s(r,t) = \frac{Q_0}{4 \pi kD} W(u), \,\,\,\,\, u = \frac{r^2 S}{4 kD t} $$

We kunnen dus de $Q(t)$ berekenen die de verlaging $s(\rho,t)$ de gewenste waarde geeft. Om de verlaging op de gewenste waarde te houden moet $Q$ op een bepaalde wijze in de tijd variëren. Hoe we deze Q(r,t) precies kunnen berekenen weten we niet, maar we kunnen de constante $Q$ die nodig is om de verlaging $s(\rho,t)$ te bereiken wel berekenenen voor elk tijdstip $t$. Deze in de tijd variërende $Q_0$ kunnen we als de gezochte benadering beschouwen.

Dit wordt hieronder uitgewerkt en grafisch geverifieerd. Het is de bedoeling om de verlaging op de rand van de bouwput, $r=\rho$ vast te houden na een zekere aanlooptijd $t_1$, die nodig wordt geacht om de gewenste verlaging dan voor het eerst te behalen.

Conclusie:

alleen wanneer de radius van de straal waarop de gewenste verlaging moet geleden groot wordt zijn de afwijkingen aanzienlijk. In de meeste gevallen is het verschil verwaarloosbaar voor de praktijk.

De aanloopperiode is die waarin de grond onder de te realiseren bouwput wordt leeggepompt met constant debiet om de bemaling in eerste instantie te realiseren. Dan volgt een periode waarin de verlaging wordt gehandhaafd, waarbij het debiet geleidelijk in de tijd daalt. De laatste tak van de grafiek is de periode na afloop van de bemaling.

Ondat de bemaling wordt berekend alsof die met een enkele put wordt uitgevoerd in het midden van de bouwput is de berekende verlaging binnen de bouwput groter dan die ter plekke van de omtrek, zijnde de gewenste verlaging en de verlaging die binnen de gehele bouwput zou worden bereikt wanneer die zou worden uitgevoerd met putten langs de omtrek van de bouwput. Dus in feite verlaagt de bemaling met een enkele central put de grondwaterstand binnend de bouwput veel meer dan de bedoeling was en is het uitgemalen debiet dus te hoog. Een bemaling langs de omtrek van de bouwput kan dus met een wat lagere onttrekking af. Op lange termijn is het volume dat binnen de bouput teveel wordt uitgemalen met de centrale put echter verwaarloosbaar en valt uiteindelijk volledig weg bij een stationaire onttrekking volgens Dupuit of Hantush.

Omdat we het onttrekkingsdebiet op elk moment gelijk stellen aan het debiet dat nodig is om pas op dat moment de gewenste verlaging te bereiken, malen we met deze berekeningsmethode wat teveel af, waardoor de bereende verlaging uiteindelijk wat hoger is dan de gwenste verlaing. Maar het verschil is voor de praktijk verwaarloosbaar.

Analytische oplossingen voor verloop van onttrekking bij constante verlaging

Er bestaan analytische formules voor het debiet tussen twee ringen wanneer het peil op de binnenste ring op $t=0$ met een vaste waarde verandert (Zie Carslaw & Jaeger (1959)). Er is ook een rapport van de USGS (Bennet & Patten (1964)) waarin dit vraagstuk wordt behandeld en analytisch wordt opgelost. Echter de formule bevat zeer slecht convergerende besselfuncties is daardoor nauwelijks te integreren. Ook Bruggeman (1999) heeft een analytische oplossing voor dit probleem, die wat beter is te hanteren, maar toch nog redelijk complex. Mijn conclusie is dat de hiervoor voorgestelde benadering veel praktischer is en bovendien nauwkeurig genoeg. Onnauwkeurigheden door onvoldoende kennis van de bodemconstanten of door variërend debiet of minder goed afgesteld debie, of door vooraf onbekende neerslag tijdens de bemaling zullen groter zijn dan de fout die met de hiervoor uitgewerkte berekeningsmethode wordt gemaakt.

Bruggeman (1999) geeft een analytische oplossing (nr. 223.02) voor de stroming buiten een cylinder met straal $R$ waar het peil op $t=0$ plotseling met een constante waarde verandert.

$$ \phi(r,t)=h\left\{1 - \frac 2 \pi \int_0^\infty \frac 1 u h(u, r) \exp\left(-\frac{u^2 t}{\beta^2 R^2}\right) du \right\} $$

terwijl

$$ h(u,r) = \frac{J_0(u)Y_0(\frac r R u) - Y_0(u)J_0(\frac r R u)}{J_0^2(u) + Y_0^2(u)} $$

We willen echter een relatie voor het debiet op afstand $r=R$.

$$ Q|_{r=R} = -2 \pi r T \frac{\partial \phi(r, t)}{\partial r}|_{r=R}$$

partial differentiëren naar $r$

$$\frac{\partial h(u,r)}{\partial r}|_{r=R} = \frac u R \frac{Y_0(u)J_1(u) - J_0(u)Y_1(u) }{J_0^2(u) + Y_0^2(u)}= \frac u R g(u)$$

met

$$g(u) = \frac{Y_0(u)J_1(u) - J_0(u)Y_1(u) }{J_0^2(u) + Y_0^2(u)}$$

zodat

$$ \frac{\partial \phi(r,t)}{\partial r} = -h\frac 2 \pi \int_0^\infty \frac 1 u \frac{\partial h(u, r)}{\partial r} \exp\left(-\frac{u^2 t}{\beta^2 R^2}\right) du = -\frac {2 h} {\pi R} \int_0^\infty g(u) \exp\left(-\frac{u^2 t}{\beta^2 R^2}\right) du $$$$Q_{r=R}(t) = -2 \pi R T \frac{\partial \phi(r, t)}{\partial t} =-2 \pi R T \,\, \left\{- \frac {2 h} {\pi R} \int_0^\infty g(u) \exp\left(-\frac{u^2 t}{\beta^2 R^2}\right) du\right\}$$

zodat tenslotte

$$Q_{r=R}(t) = 4 T \,\, h \int_0^\infty g(u) \exp\left(-\frac{u^2 t}{\beta^2 R^2}\right) du$$

Deze oplossing valt wel redelijk (numeriek) te integreren. Het is echter essentieel om hem te vergelijken met de uitkomst van een numeriek zoals Modflow of analytisch model zoals TTIM dat deze situatie wel goed kan berekenen.

Voorbeeld

Hierna wordt de oplossing van Bruggeman vergeleken met het resultaat van TTIM.

Uit de tabel hieronder blijkt dat alle afleidingen consistent en dus waarschijnlijk correct zijn, maar dat het resulterende debiet veel kleiner is dan volgens TTIM en de vereenvoudigde oplossing die hiervoor werd voorgesteld en afgeleid, terwijl die twee vlak bij elkaar liggen. Mogelijk zit er een fout in de formule die door Bruggeman wordt gegeven.

Use finite differende model to compute the flow at the face of the cylinter

With a fnite difference model it is straightforward to compute the flow at a cylinder with given radius $R$ caused by a sudden head change at this outer face of the cylinder (there is no flow inside the cylinder). We can do the same simulation using TTIM by applying a well with the same radius $R$.

Er is een klein verschil tussen de stroming die ttim berekent en die van het eindige differentiemodel. Het is me nog niet duidelijk hoe dat komt. In elk geval is de uitkomst van de formule van het boek van Bruggeman (1999), oplossing 223.02 veel kleiner en duidelijk niet goed (zie hierna).

De uitkomst volgens de formule van Bruggeman (1999) is duidelijk veel te laag en neemt te snel af in de tijd. Het is me nog niet duidelijk of er een fout zit in het boek of in mijn implementatie. Het numerieke model geeft een iets te groot debiet. Het is me nu nog niet duidelijk waar dat verschil met ttim vandaan komt. ttim is ongetwijfeld goed. Dat blijkt ook uit de eenvoudige benadering met convolutie, die lijkt zeer goed op die van ttim.

Gebruik van TTIM

TTIM is een open source API voor het analytische modelleren van grondwaterstroming in meerdere lagen. Uitgangspunt is een initiële situaite met de stijghoogten overal gelijk aan nul, die zich aanpast onder invloed van stresses die met analytische elementen worden opgelegd. Zulke analytische elementen zijn putten, "line sinks", "line doublets" en een "circular area sink", een cirkelvorming element met opgelegd neerslagoverschot (voeding). De line sinks kunnen aan elkaar worden geregen tot complexe sloot- beek en rivierconfiguraties of lijnonttrekkingen, met gegeven onttrekking per lengteeenheid of met geegven stijghoogte die al dan niet via een weerstand het grondwater in de aquifers beïnvloedt. De line doublets kunnen worden gebruikt om weerstand tegen horizontale stroming aan te brengen, zoals damwanden en gesloten randen, waarbij het grondwater aan een zijde van de dipool al dan niet via een weerstand communiceert met die aan de andere zijde. Lek tussen aquifers is automatisch besloten in de oplossing omdat het meerlagensysteem tussen elk van de lagen een weerstand krijgt opgelegd.

Voor de specifieke onttreking $\sigma_i$ [L^3/T/L] = [L^2/T] uit ene line-sink in laag $i$ geldt

$$ \sigma_i = w(h_i - h_{ls}) / c$$

met $w$ [L] de breedte van de line-sink, en $c$ [T] de weerstand ervan.

De stroming door (loodrecht op) een doublet is $$ q_n = (h^- - h^+) / c $$

Voor putten geldt dat de onttrekking uit laag $i$, $Q_i$ gelijk is aan

$$ Q_i = 2 \pi r_w H_i(h_i - h_w) / c $$

met $r_w$ de radius van de put, $c$ de intredeweerstand van de put, $h_w$ de stijghoogte in de put en $h_i$ die in laag $i$ aan de buitenzijde van de put.

Het element van type Well heeft een opgelegd debiet die van type HeadWell een opgelegde stijghoogt in de put. De onttrekking per laag wordt berekend op basis van het doorlaatvermogen van de lagen en de stijghoogte ter hoogte van elke laag.

Wat moeilijker verloopt is het ruimtelijk aanpassen van het doorlaatvermogen of de berginngscoëfficiënten van de watervoerende lagen of van de weerstand tussen de lagen. Theoretisch is dit wel mogelijk, het programma MLAEM had deze mogelijkheid bijvoorbeeld 20 jaar geleden al, maar TTIM heeft dat nog niet. TTIM is mathematisch wellicht complexer dan het oudere MLAEM omdat TTIM direct analytische meerlagensystemen implementeert, terwijl MLAEM deze systemen element per element opbouwt, en daarmee veel complexer is wat betreft de invoer.

De mogelijkheden van ttim zijn dus enigszins beperkt. Tegelijkertijd is het een zeer krachtig greedschap voor veel situaties.

Put met constant debiet (Vergelijk Theis met ttim uitkomst (example 1 in Readthedocs ttim))

Voorbeeld 2

Put met readius r met gegeven verlaging op deze radius

TTIM omgaan met gebiedsweerstand of slotenpatroon

Het beken- en slotenpatroon in een gebied is een uiting van het grondwatersysteem. Een hooggelegen gebied is normaliter een infiltratiegebied met weinig open waterlopen. Het omgekeerde geldt voor en relatie laag gelegen gebied, dat normaliter wordt gekenmerkt door kwel met veel drainagemiddelen. En dan zijn er overgangsgebieden die tussen beide in liggen. Een hoog gelegen gebied heeft hierdoor een hoge drainageweerstand en een laag gelegen gebied een kleine, terwijl die in een overgangsgebied ruimtelijk verloopt. De drainageweerstand is de gemiddlde grondwaterstand $h$ [L] ten opzichte van het peil in het oppervlaktwater $h_{ow}$ gedeeld door de gemiddelde grondwateraanvulling $N$ [L/T]

$$ c_{dr} = \frac {\overline h - h_{ow}} N $$

of

$$ N = \frac{ \overline h - h_{ow}} {c_{dr}} $$

De drainageweerstand kan uitgedrukt worden in slootbreedte $w$, de slootafstand $L$ en slootbodemweerstand $c_{sl}$ [L/T]

$$ c_{sl} = \frac w L c_{dr} $$

of

$$ c_{dr} = \frac L w c_{sl} $$

In ttim kunnen we de gewenste gebieds of drainageweerstand krijgen door de slootafstand te kiezen. Nemen we de slootbodemweerstand gelijk aan $c_{sl} = 1$ d en de slootbreedte $w = 1$ m, dan volgt $c_{dr} = L$. Let op dat de dimensie van $L$ nu dagen is geworden, wat verwarrend kan zijn. Beter is het om $w$ en $c_{sl}$ eenvoudig in de formules te laten staan, dan is verwarring onmogelijk:

$$ c_{dr} = L \frac {c_{sl}} w $$

Door parallelle sloten in het model op te nemen met een gegeven breedte en weerstand is het mogelijk om een gebiedsweerstand te simuleren. Men kan eenvoudig de slootbreedte $w$ gelijk aan 1 m nemen en de sloot weerstand gelijk aan $L \gamma$ om de gewenste gebiedsweerstand in het model te creëren en deze indien gewenst ruimtelijke te variëren.

Uiteraard kunnen waterlopen ook worden ingelezen uit een geografisch informatie systeem of een andere database en worden omgezet in de analytische elementen waar ttim mee werkt. De drainageweerstand is daar dan de resultante van. Bij een grotere dichtheid aan open waterlopen is het veelal gemakkelijker om met de wat abstracte drainageweerstand te werken. Het slotenpatroon kan dan in een model worden vervangen door een verticale gebiedsweerstand. Omgekeerd kan een verticale gebiedsweerstand wroden gesimuleerd in ttim door het plaatsen van een aantal lijnelementen op beperkte onderlinge afstand, zoals hiervoor is uitgelegd. Hiermee kunnen situaties worden geschematiseerd in ttim waarin het ene gebiedstype geleidelijk of abrupt overgaat in het andere, iets dat veelvuldig voorkomt in de praktijk.

In het voorbeeld hierna wordt een dergelijke, ruitmelijk variërende drainageweerstand gesimuleerd door middel van een aantal parallelle sloten. We gebruiken hiervoor het analytische element HeadLineSinkString, dat bestaat uit een aantal aan elkaar geregen HeadLineSinks, LineSinks dus met opgelegde stijghoogte die een bodemweerstand en een breedte hebben.

Voorbeeld

Stel de drainageweerstand ten oosten van de put, in een hydrologisch overgangsgebied verloopt van 1000 d af naar 100 d. We kunnen dit modelleren met een aantal evenwijdige HeadLineSinkStrings (sloten) in N-Z richting ten oosten van de put, met een weerstand die van sloot tot sloot verschilt. We leggen in dit voorbeeld de sloten op 200 m uit elkaar, dus $L=200$ m en passen de slootweerstand $c_{sl} = L / w$ aan zodat de gewenste in de ruimte verlopende drainageweerstand wordt gesimuleerd. We kiezen gemakshalve $w=1$ m.

ttim model

Stijghoogte langs een doorsnede, in dit geval van west naar oost door de put.

Toon resultaten in contourplot

Het berekenen van de punten voor de contourplot kost veel tijd. In dit geval ca. 1m13.5 s. Dus vergt wat geduld.

De weerstand van de sloten neemt van west naar oost af om het overgangsgebied tussen het hoge gebied links en het lage gebied rechts te modelleren.

TTIM: Same but packed in a few funtions for more easily computing scenarios

The following few cells do the same same as above, but now uses some functions to more easily compute scenarios.

Het actuele drainage-scenario en de verlagingscontouren